Tarea 9 - Modelos Lineales
Shanthal Chavarría
2021-12-21
## [1] "es_ES.UTF-8/es_ES.UTF-8/es_ES.UTF-8/C/es_ES.UTF-8/en_US.UTF-8"1. [20 puntos] El archivo covid19_CR.csv contiene la cantidad de casos nuevos diarios registrados en Costa Rica desde el 6 de marzo del 2020 hasta el 23 de abril de 2021. Con la tabla realice lo siguiente:
'data.frame': 414 obs. of 2 variables:
$ FECHA : chr "2020-03-06" "2020-03-07" "2020-03-08" "2020-03-09" ...
$ casos_nuevos: int 2 5 3 2 1 9 1 3 1 8 ...a) Convierta a serie de tiempo en un objeto tipo ts utilizando una frecuencia de forma que el Patrón-Estacional sea semanal.
library(lubridate)
library(dplyr)
datos <- datos%>%
mutate(FECHA = as.POSIXct(FECHA, format = "%Y-%m-%d"))
fecha.inicio <- datos[1,1]
fecha.final <- datos[nrow(datos),1]
serie <- ts(datos$casos_nuevos, frequency = 7, start =c(1,wday(fecha.inicio)))
serieTime Series:
Start = c(1, 6)
End = c(60, 6)
Frequency = 7
[1] 2 5 3 2 1 9 1 3 1 8 6 9 19 18 26 4 17 24 19 24 30 32 32 19 16
[26] 17 28 21 20 19 19 13 16 19 37 19 19 18 17 6 8 16 7 6 5 2 7 12 5 1
[51] 6 2 2 8 8 6 6 8 6 3 13 6 4 8 7 12 9 3 11 15 13 10 10 3 16
[ reached getOption("max.print") -- omitted 339 entries ](b) Utilizando el último mes de abril para pruebas y el resto de fechas para entrenamiento genere los modelos de HOLT-WINTERS, HOLT-WINTERS Calibrado y Redes Neuronales, luego en un solo gráfico muestre la serie de entrenamiento, la serie de prueba y el resultado de la predicción de cada uno de los modelos anteriores.
library(xts)
library(dplyr)
library(ggplot2)
library(forecast)
library(dygraphs)
library(lubridate)
train.ts <- head(serie, -23)
test.ts <- tail(serie, 23)HoltWinters
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
57.57143 559.6448 411.04690 708.2426 332.38393 786.9056
57.71429 519.7350 367.15562 672.3144 286.38497 753.0850
57.85714 525.4608 368.79858 682.1230 285.86661 765.0550
58.00000 312.7390 151.89525 473.5828 66.74970 558.7284
58.14286 265.0959 99.97441 430.2174 12.56437 517.6274
58.28571 502.9149 333.42202 672.4078 243.69790 762.1319
58.42857 554.6465 380.69088 728.6022 288.60432 820.6887
58.57143 551.4222 356.06848 746.7759 252.65446 850.1899
58.71429 511.5124 311.91110 711.1137 206.24855 816.7763
58.85714 517.2382 313.29415 721.1823 205.33267 829.1438
59.00000 304.5165 96.13645 512.8965 -14.17326 623.2062
59.14286 256.8733 43.96623 469.7804 -68.73997 582.4866
59.28571 494.6923 277.16895 712.2157 162.01902 827.3657
59.42857 546.4240 324.19689 768.6510 206.55700 886.2909
59.57143 543.1996 301.42933 784.9699 173.44386 912.9554
[ reached 'max' / getOption("max.print") -- omitted 8 rows ] Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
57.57143 559.6448 411.04690 708.2426 332.38393 786.9056
57.71429 519.7350 367.15562 672.3144 286.38497 753.0850
57.85714 525.4608 368.79858 682.1230 285.86661 765.0550
58.00000 312.7390 151.89525 473.5828 66.74970 558.7284
58.14286 265.0959 99.97441 430.2174 12.56437 517.6274
58.28571 502.9149 333.42202 672.4078 243.69790 762.1319
58.42857 554.6465 380.69088 728.6022 288.60432 820.6887
58.57143 551.4222 356.06848 746.7759 252.65446 850.1899
58.71429 511.5124 311.91110 711.1137 206.24855 816.7763
58.85714 517.2382 313.29415 721.1823 205.33267 829.1438
59.00000 304.5165 96.13645 512.8965 -14.17326 623.2062
59.14286 256.8733 43.96623 469.7804 -68.73997 582.4866
59.28571 494.6923 277.16895 712.2157 162.01902 827.3657
59.42857 546.4240 324.19689 768.6510 206.55700 886.2909
59.57143 543.1996 301.42933 784.9699 173.44386 912.9554
[ reached 'max' / getOption("max.print") -- omitted 8 rows ]calibrar.HW <- function(entrenamiento, prueba, paso = 0.1) {
# se calculan todas las combinaciones para los parametros
params <- purrr::cross(list(a = seq(0, 1, by = paso), b = seq(0, 1, by = paso), g = seq(0, 1, by = paso)))
# se calcula un modelos para cada combinacion de parametros
hw_secure <- purrr::possibly(forecast::hw, otherwise = NULL)
models <- invisible(purrr::map(params, ~suppressWarnings(hw_secure(
entrenamiento, alpha = ifelse(.$a == 0, F, .$a), beta = ifelse(.$b == 0, F, .$b),
gamma = ifelse(.$g == 0, F, .$g), h = length(prueba)))))
# se realiza la prediccion con cada uno de los modelos
predictions <- purrr::map(models, ~{
if (is.null(.)) {
return(NULL)
}
forecast(., h = length(prueba))
})
# se calcula el error para cada prediccion
error <- purrr::map_dbl(predictions, ~{
if (is.null(.)) {
return(Inf)
}
sum((as.numeric(prueba) - as.numeric(.$mean))^2)
})
# se retorna el modelo con el menor error
best_model <- models[[which.min(error)]]
p <- params[[which.min(error)]]
best_model$call <- call("HoltWinters", y = quote(datos),
alpha = ifelse(p$a == 0, F, p$a),
beta = ifelse(p$b == 0, F, p$b),
gamma = ifelse(p$g == 0, F, p$g))
return(best_model$call)
}
calibrar.HW(train.ts, test.ts)[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"
[1] "Model: ETS(A,A,A)"HoltWinters(y = datos, alpha = 0.4, beta = 0.3, gamma = 0.1)### calibrado HW
modelo.calibrado <- hw(train.ts, alpha = 0.4, beta = 0.3, gamma = 0.1, h = 23)
modelo.calibrado Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
57.57143 585.3225 381.026468 789.6185 272.8787 897.7662
57.71429 638.3513 388.976271 887.7264 256.9651 1019.7376
57.85714 669.7678 347.394187 992.1415 176.7399 1162.7958
58.00000 482.9130 65.228790 900.5971 -155.8798 1121.7058
58.14286 410.1700 -120.213309 940.5532 -400.9813 1221.3212
58.28571 650.8701 -6.379086 1308.1192 -354.3058 1656.0459
58.42857 755.0108 -41.219632 1551.2412 -462.7185 1972.7401
58.57143 785.8356 -171.308027 1742.9793 -677.9893 2249.6606
58.71429 838.8645 -276.187504 1953.9165 -866.4604 2544.1894
58.85714 870.2810 -412.073131 2152.6352 -1090.9104 2831.4724
59.00000 683.4261 -774.966899 2141.8191 -1546.9935 2913.8458
59.14286 610.6831 -1031.963445 2253.3297 -1901.5280 3122.8943
59.28571 851.3833 -983.304048 2686.0705 -1954.5288 3657.2953
59.42857 955.5239 -1078.632386 2989.6803 -2155.4497 4066.4976
59.57143 986.3488 -1263.042021 3235.7396 -2453.7976 4426.4952
[ reached 'max' / getOption("max.print") -- omitted 8 rows ] Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
57.57143 585.3225 381.026468 789.6185 272.8787 897.7662
57.71429 638.3513 388.976271 887.7264 256.9651 1019.7376
57.85714 669.7678 347.394187 992.1415 176.7399 1162.7958
58.00000 482.9130 65.228790 900.5971 -155.8798 1121.7058
58.14286 410.1700 -120.213309 940.5532 -400.9813 1221.3212
58.28571 650.8701 -6.379086 1308.1192 -354.3058 1656.0459
58.42857 755.0108 -41.219632 1551.2412 -462.7185 1972.7401
58.57143 785.8356 -171.308027 1742.9793 -677.9893 2249.6606
58.71429 838.8645 -276.187504 1953.9165 -866.4604 2544.1894
58.85714 870.2810 -412.073131 2152.6352 -1090.9104 2831.4724
59.00000 683.4261 -774.966899 2141.8191 -1546.9935 2913.8458
59.14286 610.6831 -1031.963445 2253.3297 -1901.5280 3122.8943
59.28571 851.3833 -983.304048 2686.0705 -1954.5288 3657.2953
59.42857 955.5239 -1078.632386 2989.6803 -2155.4497 4066.4976
59.57143 986.3488 -1263.042021 3235.7396 -2453.7976 4426.4952
[ reached 'max' / getOption("max.print") -- omitted 8 rows ]Redes Neuronales
Series: train.ts
Model: NNAR(22,1,40)[7]
Call: nnetar(y = train.ts, size = 40)
Average of 20 networks, each of which is
a 22-40-1 network with 961 weights
options were - linear output units
sigma^2 estimated as 97.17 Point Forecast
57.57 649.2158
57.71 712.5163
57.86 506.9182
58.00 424.1000
58.14 341.1344
58.29 474.8286
58.43 579.1668
58.57 625.3498
58.71 603.6568
58.86 535.3503
59.00 452.2606
59.14 408.6931
59.29 541.1187
59.43 619.4778
59.57 690.0131
59.71 667.1313
59.86 542.8352
60.00 511.9385
60.14 497.6412
60.29 560.1245
60.43 768.5957
60.57 763.6815
60.71 804.0039predicciones <- ts.union(
train.ts, Original = test.ts,
`HW` = pred.HW$mean,
`HW Calibrado` = pred.calibrado$mean,
`Redes neuronales` = pred.RN$mean
)
prediccionesTime Series:
Start = c(1, 6)
End = c(60, 6)
Frequency = 7
train.ts Original HW HW Calibrado Redes neuronales
1.714286 2 NA NA NA NA
1.857143 5 NA NA NA NA
2.000000 3 NA NA NA NA
2.142857 2 NA NA NA NA
2.285714 1 NA NA NA NA
2.428571 9 NA NA NA NA
2.571429 1 NA NA NA NA
2.714286 3 NA NA NA NA
2.857143 1 NA NA NA NA
3.000000 8 NA NA NA NA
3.142857 6 NA NA NA NA
3.285714 9 NA NA NA NA
3.428571 19 NA NA NA NA
3.571429 18 NA NA NA NA
3.714286 26 NA NA NA NA
[ reached getOption("max.print") -- omitted 399 rows ]c) Con un gráfico mida el error de cada uno de los modelos anteriores y determine cual de los modelos es el mejor.
RSS <- function(Pred, Real) {
return(sum((Real - Pred)^2))
}
MSE <- function(Pred, Real) {
N <- length(Real)
rss <- sum((Real - Pred)^2)
return((1/N) * rss)
}
RMSE <- function(Pred, Real) {
N <- length(Real)
rss <- sum((Real - Pred)^2)
return(sqrt((1/N) * rss))
}
RE <- function(Pred, Real) {
res <- sum(abs(Real - Pred))/sum(abs(Real))
return(res)
}
tabla.errores <- function(predicciones, real, nombres = NULL) {
r <- data.frame()
for (pred in predicciones) {
r <- rbind(r, data.frame(MSE = MSE(pred, real), RMSE = RMSE(pred, real),
RE = RE(pred, real), CORR = cor(pred, real)))
}
row.names(r) <- nombres
return(r)
}
grafico.errores <- function(errores) {
library(ggplot2)
library(reshape)
centros <- as.data.frame(apply(errores, 2, function(i) scales::rescale(i, to = c(0,
100))))
res <- melt(t(centros), varnames = c("E", "Modelos"))
res <- res[order(res$E, decreasing = F), ]
res$M <- as.character(res$M)
y = c(0, 25, 50, 75, 100)
ggplot(res, aes(x = E, y = value, group = Modelos, color = Modelos, fill = Modelos)) +
geom_polygon(alpha = 0.3, size = 1) + geom_point(size = 3) + theme_minimal() +
theme(axis.text.y = element_blank()) + xlab("") + ylab("") + scale_y_continuous(limits = c(-10,
100), breaks = y) + annotate("text", x = 0.5, y = y, label = paste0(y, "%"),
color = "gray60") + ggproto("CordRadar", CoordPolar, theta = "x", r = "y",
start = 0, direction = sign(1))
}errores <- tabla.errores(predicciones = list(pred.HW$mean, pred.calibrado$mean,pred.RN$mean), real = test.ts)
row.names(errores) <- c("HW", "HW.Calibrado","Redes Neuronales")
errores
El mejor modelo es Holt -Winters Calibrado, debido a que tiene mayor correlación y el menor error de todos.
d) Con el mejor modelo encontrado en el punto anterior genere la predicción de un mes, pero esta vez utilizando toda la serie de tiempo. Grafique la serie original y la predicción junto con el limite inferior y superior.
#con toda la serie
modelo <- hw(serie, alpha = 0.4, beta = 0.3, gamma = 0.1, h = 30, level = c(95))
modelo Point Forecast Lo 95 Hi 95
60.85714 1989.630 1658.9488 2320.312
61.00000 2047.884 1644.2364 2451.533
61.14286 2258.683 1736.8770 2780.490
61.28571 2825.057 2148.9767 3501.136
61.42857 3169.665 2311.1663 4028.164
61.57143 3385.451 2321.6016 4449.300
61.71429 3558.011 2269.2021 4846.821
61.85714 3845.539 2296.2690 5394.809
62.00000 3903.793 2098.9267 5708.659
62.14286 4114.592 2038.9239 6190.260
62.28571 4680.965 2320.3536 7041.576
62.42857 5025.574 2366.7224 7684.425
62.57143 5241.359 2271.6632 8211.055
62.71429 5413.920 2121.3558 8706.484
62.85714 5701.447 2060.4963 9342.398
63.00000 5759.702 1774.4760 9744.927
63.14286 5970.501 1630.1258 10310.875
63.28571 6536.874 1830.8459 11242.901
63.42857 6881.482 1799.6286 11963.336
63.57143 7097.268 1629.7126 12564.823
63.71429 7269.829 1406.9673 13132.690
63.85714 7557.356 1278.0649 13836.647
64.00000 7615.610 923.2521 14307.968
64.14286 7826.409 711.9922 14940.826
64.28571 8392.782 847.5190 15938.045
[ reached 'max' / getOption("max.print") -- omitted 5 rows ]pred <- forecast(modelo, h = 30)
p <- ts.union(
prediccion = pred$mean, LimInf = pred$lower, LimSup = pred$upper
)
fecha <- fecha_inicio + days(0:443)
predicciones <- ts.union(serie, p)
prediccionesTime Series:
Start = c(1, 6)
End = c(65, 1)
Frequency = 7
serie p.prediccion p.LimInf p.LimSup
1.714286 2 NA NA NA
1.857143 5 NA NA NA
2.000000 3 NA NA NA
2.142857 2 NA NA NA
2.285714 1 NA NA NA
2.428571 9 NA NA NA
2.571429 1 NA NA NA
2.714286 3 NA NA NA
2.857143 1 NA NA NA
3.000000 8 NA NA NA
3.142857 6 NA NA NA
3.285714 9 NA NA NA
3.428571 19 NA NA NA
3.571429 18 NA NA NA
3.714286 26 NA NA NA
3.857143 4 NA NA NA
4.000000 17 NA NA NA
4.142857 24 NA NA NA
[ reached getOption("max.print") -- omitted 426 rows ]predicciones <- xts(xts(predicciones, order.by = fecha))
dygraph(predicciones, width = "100%") %>%
dySeries("serie", label = "Actual") %>%
dySeries(c("p.LimInf", "p.prediccion", "p.LimSup"), label = "Predicción") %>%
dyRangeSelector(height = 20, strokeColor = "") %>%
dyOptions(axisLineColor = "navy", gridLineColor = "lightblue")2. [35 puntos] La tabla trafico trafico_tren_V2.csv contiene la cantidad de personas que se suben al metro por hora. Tomando los datos a partir del 15 de Enero del 2016 realice lo siguiente:
datos <- read.table("../datos/trafico_tren_V2.csv", sep = ",", dec = ".", header = T)
fecha.inicio <- as.POSIXct("2016-01-15 00:00:00")
datos <- datos %>%
mutate(date_time = as.POSIXct(date_time, format = "%Y-%m-%d %H:%M:%S"))%>%
filter(date_time >= fecha.inicio)
fecha.final <- datos[nrow(datos),1](a) Encuentre las fechas faltantes y con un suavizado de 25 impute valor a dichas fechas.
total.fechas <- seq(fecha.inicio, fecha.final, by = "hour")
faltan.fechas <- total.fechas[!total.fechas %in% datos$date_time]
datos <- union_all(datos, data.frame(date_time = faltan.fechas))
datos <- datos[order(datos$date_time), ]
# función suavizado
suavizado <- function(datos, n) {
library(zoo)
if(n %% 2 == 0) {
izquierda = rep(NA, (floor(n/2) - 1))
derecha = rep(NA, floor(n/2))
} else {
izquierda = derecha = rep(NA, floor(n/2))
}
datos <- c(izquierda, datos, derecha)
return(rollapply(datos, n, mean, na.rm = T))
}
trafico.suavizado <- suavizado(datos$traffic_volume,25)
datos$traffic_volume[which(is.na(datos$traffic_volume))] <- trafico.suavizado[which(is.na(datos$traffic_volume))]
datos(b) Convierta a serie de tiempo (utilice Patrón-estacional diario).
Time Series:
Start = c(1, 1)
End = c(990, 24)
Frequency = 24
[1] 753.000 481.000 3515.375 3515.375 829.000 3845.111 5394.000 3940.500
[9] 5502.000 3849.727 4447.000 3643.923 5154.000 3508.429 5563.000 3526.786
[17] 6483.000 3516.714 4799.000 3128.067 2942.000 2934.800 3324.000 1700.000
[25] 1467.000 1033.000 738.000 2877.750 2981.176 688.000 1211.000 1870.000
[33] 2603.000 2691.000 3746.000 4207.000 4601.000 2628.812 4632.000 2579.812
[41] 4636.000 2455.471 3917.000 2482.529 2683.000 2484.824 2349.000 2082.000
[49] 1065.000 2333.588 644.000 343.000 466.000 2229.294 929.000 1430.000
[57] 1988.000 2524.000 2857.000 3172.000 2006.059 3954.000 3746.000 1940.706
[65] 3749.000 2144.188 2891.000 2457.059 2084.000 2679.882 1284.000 977.000
[73] 2816.062 335.000 263.000
[ reached getOption("max.print") -- omitted 23685 entries ](c) Usando las últimas 24 horas para pruebas y el resto de fechas para entrenamiento genere los modelos de HOLT-WINTERS, HOLT-WINTERS Calibrado y Redes Neuronales, luego en un solo gráfico muestre la serie de entrenamiento, la serie de prueba y el resultado de la predicción de cada uno de los modelos anteriores.
HOLT-Winter
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
990.0000 2793.276 1922.0662 3664.486 1460.87565 4125.676
990.0417 2474.484 1357.4468 3591.521 766.12303 4182.845
990.0833 2432.369 1114.5446 3750.193 416.93056 4447.807
990.1250 2357.611 865.7421 3849.479 75.99452 4639.227
990.1667 2714.851 1067.1839 4362.518 194.96162 5234.740
990.2083 4202.956 2412.9655 5992.945 1465.40210 6940.509
990.2500 6156.367 4234.5330 8078.200 3217.17573 9095.558
990.2917 6760.799 4715.5735 8806.025 3632.89647 9888.702
990.3333 6560.519 4398.9052 8722.133 3254.61593 9866.422
990.3750 6328.135 4056.0602 8600.210 2853.29647 9802.974
990.4167 6117.080 3739.6447 8494.514 2481.10674 9753.052
990.4583 6352.817 3874.4720 8831.161 2562.51564 10143.118
990.5000 6557.475 3982.1477 9132.802 2618.85178 10496.098
990.5417 6532.777 3863.9659 9201.588 2451.18263 10614.371
990.5833 6720.864 3961.7129 9480.016 2501.10634 10940.623
[ reached 'max' / getOption("max.print") -- omitted 9 rows ] Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
990.0000 2793.276 1922.0662 3664.486 1460.87565 4125.676
990.0417 2474.484 1357.4468 3591.521 766.12303 4182.845
990.0833 2432.369 1114.5446 3750.193 416.93056 4447.807
990.1250 2357.611 865.7421 3849.479 75.99452 4639.227
990.1667 2714.851 1067.1839 4362.518 194.96162 5234.740
990.2083 4202.956 2412.9655 5992.945 1465.40210 6940.509
990.2500 6156.367 4234.5330 8078.200 3217.17573 9095.558
990.2917 6760.799 4715.5735 8806.025 3632.89647 9888.702
990.3333 6560.519 4398.9052 8722.133 3254.61593 9866.422
990.3750 6328.135 4056.0602 8600.210 2853.29647 9802.974
990.4167 6117.080 3739.6447 8494.514 2481.10674 9753.052
990.4583 6352.817 3874.4720 8831.161 2562.51564 10143.118
990.5000 6557.475 3982.1477 9132.802 2618.85178 10496.098
990.5417 6532.777 3863.9659 9201.588 2451.18263 10614.371
990.5833 6720.864 3961.7129 9480.016 2501.10634 10940.623
[ reached 'max' / getOption("max.print") -- omitted 9 rows ]calibrado HW
library(xts)
library(dplyr)
library(ggplot2)
library(forecast)
library(dygraphs)
library(lubridate)
calibrar.HW <- function(serie.aprendizaje, serie.testing, paso = 0.1) {
error.c <- Inf
alpha.i <- paso # alpha no puede ser cero
while (alpha.i <= 1) {
beta.i <- 0
while (beta.i <= 1) {
gamma.i <- 0
while (gamma.i <= 1) {
mod.i <- HoltWinters(serie.aprendizaje, alpha = alpha.i, beta = beta.i,
gamma = gamma.i)
res.i <- predict(mod.i, n.ahead = length(serie.testing))
error.i <- sqrt(MSE(res.i, serie.testing))
if (error.i < error.c) {
error.c <- error.i
mod.c <- mod.i
}
gamma.i <- gamma.i + paso
}
beta.i <- beta.i + paso
}
alpha.i <- alpha.i + paso
}
return(mod.c)
}
calibrar.HW(train.ts, test.ts)Holt-Winters exponential smoothing with trend and additive seasonal component.
Call:
HoltWinters(x = serie.aprendizaje, alpha = alpha.i, beta = beta.i, gamma = gamma.i)
Smoothing parameters:
alpha: 0.1
beta : 0
gamma: 1
Coefficients:
[,1]
a -302752.86271
b -47.29718
s1 303284.38627
s2 302434.46961
s3 302123.59654
s4 302010.45111
s5 302167.58107
s6 302727.16609
s7 303778.80105
s8 304871.03307
s9 306120.50073
s10 306803.28144
s11 307110.52675
s12 307402.16449
s13 307648.60570
s14 307509.70405
s15 307565.80321
s16 307576.26255
s17 307827.64542
s18 307983.63477
s19 307653.39360
s20 306869.91626
s21 306330.01066
s22 305878.75462
s23 305475.17219
s24 306608.86271modelo.calibrado <- HoltWinters(train.ts, alpha = 0.1, beta = 0, gamma = 1)
pred.calibrado <- forecast(modelo.calibrado, h = 24)
pred.calibrado Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
990.0000 484.2264 -890.1564 1858.6092 -1617.7107 2586.163
990.0417 -412.9875 -1794.2251 968.2502 -2525.4081 1699.433
990.0833 -771.1577 -2159.2163 616.9009 -2894.0101 1351.695
990.1250 -931.6003 -2326.4465 463.2459 -3064.8335 1201.633
990.1667 -821.7675 -2223.3685 579.8334 -2965.3312 1321.796
990.2083 -309.4797 -1717.8030 1098.8436 -2463.3243 1844.365
990.2500 694.8581 -720.1556 2109.8718 -1469.2186 2858.935
990.2917 1739.7929 318.1203 3161.4655 -434.4677 3914.054
990.3333 2941.9634 1513.6629 4370.2639 757.5663 5126.361
990.3750 3577.4469 2142.5492 5012.3447 1382.9602 5771.934
990.4167 3837.3951 2395.9302 5278.8599 1632.8648 6041.925
990.4583 4081.7356 2633.7335 5529.7378 1867.2075 6296.264
990.5000 4280.8797 2826.3696 5735.3897 2056.3985 6505.361
990.5417 4094.6808 2633.6919 5555.6698 1860.2910 6329.071
990.5833 4103.4828 2636.0435 5570.9221 1859.2281 6347.737
[ reached 'max' / getOption("max.print") -- omitted 9 rows ]Redes Neuronales
Series: train.ts
Model: NNAR(43,1,10)[24]
Call: nnetar(y = train.ts, size = 10)
Average of 20 networks, each of which is
a 43-10-1 network with 451 weights
options were - linear output units
sigma^2 estimated as 189260 Point Forecast
990.00 2055.2810
990.04 1463.2901
990.08 852.0951
990.12 609.3308
990.17 427.2401
990.21 414.5691
990.25 751.5265
990.29 1165.6261
990.33 1881.2204
990.38 2499.5063
990.42 3207.6977
990.46 3799.0823
990.50 4044.5884
990.54 4134.4537
990.58 4334.4183
990.62 4308.7862
990.67 4277.5442
990.71 4285.0110
990.75 4104.9075
990.79 3768.1915
990.83 3194.2490
990.88 2778.3316
990.92 2279.5560
990.96 2071.2922predicciones <- ts.union(
train.ts, Original = test.ts,
`HW` = pred.HW$mean,
`HW Calibrado` = pred.calibrado$mean,
`Redes neuronales` = pred.RN$mean
)
prediccionesTime Series:
Start = c(1, 1)
End = c(990, 24)
Frequency = 24
train.ts Original HW HW Calibrado Redes neuronales
1.000000 753.000 NA NA NA NA
1.041667 481.000 NA NA NA NA
1.083333 3515.375 NA NA NA NA
1.125000 3515.375 NA NA NA NA
1.166667 829.000 NA NA NA NA
1.208333 3845.111 NA NA NA NA
1.250000 5394.000 NA NA NA NA
1.291667 3940.500 NA NA NA NA
1.333333 5502.000 NA NA NA NA
1.375000 3849.727 NA NA NA NA
1.416667 4447.000 NA NA NA NA
1.458333 3643.923 NA NA NA NA
1.500000 5154.000 NA NA NA NA
1.541667 3508.429 NA NA NA NA
1.583333 5563.000 NA NA NA NA
[ reached getOption("max.print") -- omitted 23745 rows ](d) Con un gráfico mida el error y determine cuál modelo es el mejor.
errores <- tabla.errores(predicciones = list(pred.HW$mean, pred.calibrado$mean,pred.RN$mean), real = test.ts)
row.names(errores) <- c("HW", "HW.Calibrado","Redes Neuronales")
grafico.errores(errores)
(e) Con el mejor modelo encontrado en el punto anterior genere la predicción de un día, pero esta vez utilizando toda la serie de tiempo. Grafique la serie original y la predicción junto con el límite inferior y superior con dygraph.
#con toda la serie
modelo <- nnetar(serie, size = 10,level = c(95))
pred <- forecast(modelo, h = 24,PI = TRUE)
p <- ts.union(
prediccion = pred$mean, LimInf = pred$lower[,
2], LimSup = pred$upper[, 2])
predicciones <- ts.union(serie, p)
horas <- fecha.inicio + hours(0:23783)
predicciones <- xts(xts(predicciones, order.by = horas))
dygraph(predicciones, width = "100%") %>%
dySeries("serie", label = "Actual") %>%
dySeries(c("p.LimInf", "p.prediccion", "p.LimSup"), label = "Predicción") %>%
dyRangeSelector(height = 20, strokeColor = "") %>%
dyOptions(axisLineColor = "navy", gridLineColor = "lightblue")3. [30 puntos] Para la tabla Cajero.csv visto en clase, genere las reglas para la predicción del retiro de dinero de los días 14, 15 y 16 de agosto del 2012. Puede usar la regla que se implementó en el ejemplo de la clase para el día 15 de agosto del 2012.
14 de agosto
cajero <- read.csv("../datos/Cajero.csv", header = T, dec = ".", sep = ";")
fechas <- seq(as.Date("1998-01-01"), as.Date("2012-07-30"), "day")
cajero$fechas <- fechas
# Eliminamos días bisiestos
cajero <- cajero[!(month(cajero$fechas) == 2 & day(cajero$fechas) == 29), ]
# verificamos si hay valores con NA
cajero$monto[is.na(cajero$monto)]numeric(0)Time Series:
Start = c(1998, 1)
End = c(1998, 6)
Frequency = 365
[1] 628000 1162000 3891000 3005000 2716000 2621000modelo.14 <- HoltWinters(serie.regla.14)
prediccion.14.agosto <- forecast(modelo.14, h = 1)
prediccion.14.agosto$mean[1] fit
11572558 [1] 16657000[1] -5084442if (error < 0) {
factor.ajuste <- as.numeric(1 + (abs(error)/prediccion.14.agosto$mean[1]))
} else {
factor.ajuste <- as.numeric(1 - (abs(error)/prediccion.14.agosto$mean[1]))
}
factor.ajuste[1] 1.439353 fit
16657000 [1] 16657000Time Series:
Start = c(1998, 1)
End = c(2012, 211)
Frequency = 365
[1] 628000 1162000 3891000 3005000 2716000 2621000 3293000 2794000
[9] 866000 2701000 2574000 2627000 3962000 11005000 1759000 2398000
[17] 5250000 3831000 3505000 3700000 3968000 2184000 1691000 2542000
[25] 2906000 2573000 4000000 6721000 4647000 2247000 7414000 5029000
[33] 3818000 3909000 4001000 3808000 1819000 3931000 2391000 2393000
[41] 2655000 4079000 3170000 1685000 4169000 7738000 6172000 3814000
[49] 3213000 2923000 985000 1567000 2403000 2737000 2757000 4695000
[57] 2584000 1376000 4354000 10641000 5647000 4938000 3869000 2513000
[65] 1074000 3463000 2389000 2616000 3201000 4733000 2194000 1153000
[73] 4310000 4764000 7648000
[ reached getOption("max.print") -- omitted 5246 entries ][1] 14539697[1] 20927760if (pred$mean[15] > max(serie)) {
pred$mean[15] <- max(serie)
}
if (pred$mean[15] < min(serie)) {
pred$mean[15] <- min(serie)
}
serie.pred <- ts.union(serie, Predicción = pred$mean)
fecha.fin <- cajero$fechas[nrow(cajero)]
# Generamos las fechas de la predicción.
fechas <- fecha.fin + days(1:30)
# Unimos con las fechas de la serie original
total.fechas <- c(cajero$fechas, fechas)
serie.pred <- xts(xts(serie.pred, order.by = total.fechas))
dygraph(serie.pred, "Pronóstico de cantidad de dinero del cajero", width = "100%") %>%
dyAnnotation("2012-08-14", text = "14") %>% # Denotamos el valor de la regla con R
dyOptions(drawPoints = TRUE, pointSize = 2) %>%
dyRangeSelector(height = 20, strokeColor = "") %>%
dyRangeSelector(fillColor = "black", strokeColor = "orange" )15 de agosto
Time Series:
Start = c(1998, 1)
End = c(1998, 6)
Frequency = 365
[1] 628000 1162000 3891000 3005000 2716000 2621000modelo.15 <- HoltWinters(serie.regla.15)
prediccion.15.agosto <- forecast(modelo.15, h = 1)
prediccion.15.agosto$mean[1] fit
13543013 [1] 15625000[1] -2081987if (error < 0) {
factor.ajuste <- as.numeric(1 + (abs(error)/prediccion.15.agosto$mean[1]))
} else {
factor.ajuste <- as.numeric(1 - (abs(error)/prediccion.15.agosto$mean[1]))
}
factor.ajuste[1] 1.153731 fit
15625000 [1] 15625000Time Series:
Start = c(1998, 1)
End = c(2012, 211)
Frequency = 365
[1] 628000 1162000 3891000 3005000 2716000 2621000 3293000 2794000
[9] 866000 2701000 2574000 2627000 3962000 11005000 1759000 2398000
[17] 5250000 3831000 3505000 3700000 3968000 2184000 1691000 2542000
[25] 2906000 2573000 4000000 6721000 4647000 2247000 7414000 5029000
[33] 3818000 3909000 4001000 3808000 1819000 3931000 2391000 2393000
[41] 2655000 4079000 3170000 1685000 4169000 7738000 6172000 3814000
[49] 3213000 2923000 985000 1567000 2403000 2737000 2757000 4695000
[57] 2584000 1376000 4354000 10641000 5647000 4938000 3869000 2513000
[65] 1074000 3463000 2389000 2616000 3201000 4733000 2194000 1153000
[73] 4310000 4764000 7648000
[ reached getOption("max.print") -- omitted 5246 entries ][1] 15216130[1] 17555328if (pred$mean[16] > max(serie)) {
pred$mean[16] <- max(serie)
}
if (pred$mean[16] < min(serie)) {
pred$mean[16] <- min(serie)
}
serie.pred <- ts.union(serie, Predicción = pred$mean)
fecha.fin <- cajero$fechas[nrow(cajero)]
# Generamos las fechas de la predicción.
fechas <- fecha.fin + days(1:30)
# Unimos con las fechas de la serie original
total.fechas <- c(cajero$fechas, fechas)
serie.pred <- xts(xts(serie.pred, order.by = total.fechas))
dygraph(serie.pred, "Pronóstico de cantidad de dinero del cajero", width = "100%") %>%
dyAnnotation("2012-08-15", text = "15") %>% # Denotamos el valor de la regla con R
dyOptions(drawPoints = TRUE, pointSize = 2) %>%
dyRangeSelector(height = 20, strokeColor = "") %>%
dyRangeSelector(fillColor = "black", strokeColor = "orange" )16 de agosto
Time Series:
Start = c(1998, 1)
End = c(1998, 6)
Frequency = 365
[1] 628000 1162000 3891000 3005000 2716000 2621000modelo.16 <- HoltWinters(serie.regla.16)
prediccion.16.agosto <- forecast(modelo.16, h = 1)
prediccion.16.agosto$mean[1] fit
13062527 [1] 9269000[1] 3793527if (error < 0) {
factor.ajuste <- as.numeric(1 + (abs(error)/prediccion.16.agosto$mean[1]))
} else {
factor.ajuste <- as.numeric(1 - (abs(error)/prediccion.16.agosto$mean[1]))
}
factor.ajuste[1] 0.709587 fit
9269000 [1] 9269000Time Series:
Start = c(1998, 1)
End = c(2012, 211)
Frequency = 365
[1] 628000 1162000 3891000 3005000 2716000 2621000 3293000 2794000
[9] 866000 2701000 2574000 2627000 3962000 11005000 1759000 2398000
[17] 5250000 3831000 3505000 3700000 3968000 2184000 1691000 2542000
[25] 2906000 2573000 4000000 6721000 4647000 2247000 7414000 5029000
[33] 3818000 3909000 4001000 3808000 1819000 3931000 2391000 2393000
[41] 2655000 4079000 3170000 1685000 4169000 7738000 6172000 3814000
[49] 3213000 2923000 985000 1567000 2403000 2737000 2757000 4695000
[57] 2584000 1376000 4354000 10641000 5647000 4938000 3869000 2513000
[65] 1074000 3463000 2389000 2616000 3201000 4733000 2194000 1153000
[73] 4310000 4764000 7648000
[ reached getOption("max.print") -- omitted 5246 entries ][1] 12642804[1] 8971169if (pred$mean[17] > max(serie)) {
pred$mean[17] <- max(serie)
}
if (pred$mean[17] < min(serie)) {
pred$mean[17] <- min(serie)
}
serie.pred <- ts.union(serie, Predicción = pred$mean)
fecha.fin <- cajero$fechas[nrow(cajero)]
# Generamos las fechas de la predicción.
fechas <- fecha.fin + days(1:30)
# Unimos con las fechas de la serie original
total.fechas <- c(cajero$fechas, fechas)
serie.pred <- xts(xts(serie.pred, order.by = total.fechas))
dygraph(serie.pred, "Pronóstico de cantidad de dinero del cajero", width = "100%") %>%
dyAnnotation("2012-08-16", text = "16") %>% # Denotamos el valor de la regla con R
dyOptions(drawPoints = TRUE, pointSize = 2) %>%
dyRangeSelector(height = 20, strokeColor = "") %>%
dyRangeSelector(fillColor = "black", strokeColor = "orange" )4. [20 puntos] Usando el Método Multiplicativo de Holt-Winters que se presenta en las siguientes ecuaciones:
\[\begin{align*} \hat{x}_{t+h|t} &= (l_t + hb_t ) \cdot s_{t-m+h^{+}_{m}} \\ l_t &= \alpha \dfrac{y_t}{s_{t-m}} + (1-\alpha)(l_{t-1} + b_{t-1}) \\ b_t &= \beta (l_t -l_{t-1}) + (1-\beta)b_{t-1} \\ s_t &= \gamma \dfrac{y_t}{(l_{t-1}+b_{t-1})}+(1-\gamma)s_{t-m} \end{align*}\]
Programe una función recursiva que permita verificar los datos que se muestran en la siguiente tabla, para esto utilice los datos de visitas de turistas a Australia. Además tome \(\alpha = 0,441\), \(\beta = 0,030\) y \(\gamma = 0,002\) (recuerde que estos parámetros se obtienen minimizado el RMSE), además tome \(l_0 = 32,49\), \(b_0 = 0,70\) y \(s_0 = 1,02\):
# Valores iniciales
m <- 4
alpha <- 0.441
beta <- 0.030
gamma <- 0.002
l0 <- 32.49
b0 <- 0.70
s0 <- 1.02
st.neg <- c(1.24,0.77,0.96,1.02)
yth <- function(t,m=4){
if(t == 1){
return((l0 + 1*b0)*st.neg[t-1 + (trunc((1-1)/m) + 1)])
}
if(t<5){
yt <- (lt(t-1) + 1*bt(t-1))*st.neg[t-1 + (trunc((1-1)/m) + 1)]
}
else{
yt <- (lt(t-1) + 1*bt(t-1))* st(t-1 - m + (trunc((1-1)/m) + 1))
}
return(yt)
}
lt <- function(t,m=4){
if(t==1){
lt <- alpha*(datos[t]/st.neg[t])+(1-alpha)*(l0 + b0)
}
if(t<5){
lt <- alpha*(datos[t]/st.neg[t])+(1-alpha)*(lt(t-1) + bt(t-1))
}
else{
lt <- alpha*(datos[t]/st(t-m))+(1-alpha)*(lt(t-1) + bt(t-1))
}
return(lt)
}
bt <-function(t){
if(t==1){
bt <- beta*(lt(1)-l0) + (1-beta)*b0
}
else{
bt <- beta*(lt(t)-lt(t-1)) + (1-beta)*bt(t-1)
}
return(bt)
}
st <- function(t, m=4){
if(t==1){
return(gamma*(datos[t] - l0 - b0) + (1-gamma)*st.neg(t))
}
if(t < 4){
st <- gamma*(datos[t] - lt(t-1) - bt(t-1)) + (1-gamma)*st.neg(t)
}
else{
st <- gamma*(datos[t] - lt(t-1) - bt(t-1)) + (1-gamma)*st.neg(t-m)
}
return(st)
}
yth(t =1)5. [20 puntos] (Suavizado exponencial simple) Se tienen \(x_1 , ... , x_n\) y se desea pronosticar (predecir) \(x_{n+h}\) con \(h \in \mathbb{N}\) y donde \(\hat{x}_{n,h}\) denota la predicción de \(x_{n+h}\). Para \(\alpha \in ]0;1[\) , se define la
predicción por suavizado exponencial simple como sigue:
\[ \hat{x}_{n,h} = \alpha \sum_{j=0}^{n-1}(1-\alpha)^{j} x_{n-j} \] Lo anterior se puede calcular recursivamente usando la siguiente fórmula:
\[ \hat{x}_{n,h} \alpha x_n + (1-\alpha)\hat{x}_{n-1,h} \] Pruebe que la predicción \(\hat{x}_{n,h}\) es asintóticamente la solución de:
\[ \hat{x}_{n,h} = arg_{x} mín \sum_{j=0}^{n-1} (1-\alpha)^j (x_{n-j} - x)^2 \]
Solución:
